capture log close
clear all
set maxvar 10000
set more off
pause off
# delimit;

local crime "";
local output "";


/* 
RD estimates (main)
*/


local exp "prwora";

local fs_vars  "age18_yes n_age18unfav n_onssi1922 n_onssi1922_exact n_onssi20 " ;
local hh_vars "n_earnhh9702 n_earnhh9712 n_emphhpost9702 n_emphhpost9712 n_inchh9702 n_inchh9712" ;
local covariates "n_male n_firstage n_fam_singmom n_fam_nopar n_dibmdr_mie n_dibmdr_mine n_dibmdr_miss n_diag1_mental n_diag1_nerv n_diag1_sense n_diag1_infec n_diag1_endo n_diag1_none n_diag1_cong n_diag1_musc n_diag1_resp n_diag1_blood n_diag1_neo n_earnhhpre9095" ;
local new_vars "*req1* *req2*" ;
local robust_vars "n_crime2char_1997_2017 n_crime2chart_1997_2017 n_crime2genr_1997_2017 n_crime2genrt_1997_2017  n_crime2nogn_1997_2017 n_crime2nognt_1997_2017 n_crime2ahinc_1997_2017" ;

local controls_all "n_male i.n_firstage i.n_diag1cat n_dibmdr_mie n_dibmdr_mine n_dibmdr_miss n_fam_singmom n_fam_nopar i.n_earnhhbins i.state" ;
local controls_no  "" ; 

local bwidth_list "264" ;
local cov_list   "no all" ; 
local demo_main "all male female" ;
local demo_oth "firstlt6 first7to13 firstgt13 remlow remhigh singno twopar parhigh parlow othmen nonment intel phys nonphys white black nonw 
				ue97high ue97low uegrhigh uegrlow ueallhigh uealllow" ;

foreach demo in `demo_main' { ;
	local vars_`demo' "n_crime* c_crime* n_1crime* n_2crime* n_3crime* n_4crime* n_yr* avg_*" ; /*  */
	local othvars_`demo' "`covariates' `fs_vars' `hh_vars'" ;
} ;
foreach demo in `demo_oth' { ;
	local vars_`demo' "n_crime* avg_*" ;
	local othvars_`demo' "`fs_vars'" ;
} ;

local poly_list "1(1)1" ;

use "`crime'\prep2\crime_prwora_prepforreg.dta", clear ;

/* Create predicted removal probability */
xi: reg n_age18unfav 
	`controls_all' i.dibdig
	if n_post == 1 ;
predict hat_removal ;

summ hat_removal, det ;
gen byte demo_q25 		= (hat_removal <  r(p25)) if hat_removal != . ;
gen byte demo_q2575 	= (hat_removal >= r(p25) & hat_removal < r(p75)) if hat_removal != . ;
gen byte demo_q75 		= (hat_removal >= r(p75) & hat_removal != .) if hat_removal != . ;
gen byte demo_remlow 	= (hat_removal <  r(p50)) if hat_removal != . ;
gen byte demo_remhigh 	= (hat_removal >= r(p50) & hat_removal != .) if hat_removal != . ;

keep if abs(dobdd_run) <= `bwidth_list' ;

/* Subsamples */
gen demo_all = 1 ;
gen byte demo_othment	= (n_diag1_mental == 1 & (substr(dibdig,1,3) != "317" & substr(dibdig,1,3) != "318" & substr(dibdig,1,3) != "319" )) ;
gen byte demo_nonment 	= (demo_othment == 0) ;
gen byte demo_intel 	= (n_diag1_mental == 1 & (substr(dibdig,1,3) == "317" | substr(dibdig,1,3) == "318" | substr(dibdig,1,3) == "319" )) ;
gen byte demo_phys 		= (n_diag1_mental != 1) ;
gen byte demo_nonphys 	= (n_diag1_mental == 1) ;
gen byte demo_twopar 	= (n_fam_singmom != 1 & n_fam_nopar != 1) ;
gen byte demo_singno 	= (n_fam_singmom == 1) | (n_fam_nopar ==1) ;
gen byte demo_nopar  	= (n_fam_nopar == 1) ;
gen byte demo_male 	 	= (n_male ==1) ;
gen byte demo_female 	= (n_male !=1) ;
gen byte demo_white		= (race == "W") ;
gen byte demo_black 	= (race == "B") ;
gen byte demo_nonw		= (race != "W") ;

egen cnty_ue97 = rowmean(cnty_urate_1997-cnty_urate_1998) ;
egen cnty_uegr = rowmean(cnty_urate_2007-cnty_urate_2012) ;
egen cnty_ueall = rowmean(cnty_urate_1997-cnty_urate_2017) ;
rename cnty_crimerate cnty_crime ;
rename cnty_share_blue cnty_dem ;
rename cnty_pop2000 cnty_pop ;

foreach cnty in crime dem pop ue97 uegr ueall { ;
	summ cnty_`cnty', det ;
	gen byte demo_`cnty'high = (cnty_`cnty' > r(p50) & cnty_`cnty' != .) ;
	gen byte demo_`cnty'low = (cnty_`cnty' <= r(p50) & cnty_`cnty' != .) ;
	replace demo_`cnty'high = . if cnty_`cnty' == . ;
	replace demo_`cnty'low = . if cnty_`cnty' == . ;
	tab demo_`cnty'high, m ;
	tab demo_`cnty'low, m ;
} ;

summ n_earnhhpre9095, det ;
replace n_earnhhpre9095 = 0 if n_earnhhpre9095 == . ;
gen byte demo_parhigh = (n_earnhhpre9095 > r(p50) & n_earnhhpre9095 != .) ;
gen byte demo_parlow = (n_earnhhpre9095 <= r(p50) & n_earnhhpre9095 != .) ;
replace demo_parhigh = . if n_earnhhpre9095 == . ;
replace demo_parlow = . if n_earnhhpre9095 == . ;
tab demo_parhigh, m ;
tab demo_parlow, m ;

gen byte demo_firstlt6 		= (n_firstage <=6) ;
gen byte demo_first7to13 	= (n_firstage > 6 & n_firstage <=13) ;
gen byte demo_firstgt13 	= (n_firstage > 13 & n_firstage < .) ;

/* RD variables */
gen dobdd_run2=dobdd_run*dobdd_run;
gen dobdd_run3=dobdd_run2*dobdd_run;
gen dobdd_run4=dobdd_run3*dobdd_run;

gen postXdobdd=n_post*dobdd_run;
gen postXdobdd2=n_post*dobdd_run2;
gen postXdobdd3=n_post*dobdd_run3;
gen postXdobdd4=n_post*dobdd_run4;

local poly_1 "dobdd_run postXdobdd";
local poly_2 "dobdd_run dobdd_run2 postXdobdd postXdobdd2";
local poly_3 "dobdd_run dobdd_run2 dobdd_run3 postXdobdd postXdobdd2 postXdobdd3";
local poly_4 "dobdd_run dobdd_run2 dobdd_run3 dobdd_run4 postXdobdd postXdobdd2 postXdobdd3 postXdobdd4";


/* ****************************************** */
/* ***** MAIN RD ESTIMATES (USE OUTREG) ***** */
/* ****************************************** */

foreach demo in all { ;

foreach bwidth in `bwidth_list' {;

	forval poly = `poly_list' {;
		
		foreach cov in `cov_list' { ;
		
				
				/* RD regressions */
				
				reg n_post n_post `poly_`poly'' 
					if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1 
					, robust;
					outreg2 using "`output'\\`exp'_bw`bwidth'_crime_poly`poly'`cov'cov_rd_`demo'.xls", e(N r2 rmse) replace;
				
				foreach var of varlist 
					`vars_`demo''
					{;
					xi: reg `var' n_post `poly_`poly''
						`controls_`cov''
						if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1
						, robust; 
					outreg2 using "`output'\\`exp'_bw`bwidth'_crime_poly`poly'`cov'cov_rd_`demo'.xls", e(N r2 rmse) append;
					} ;	
				
				foreach var of varlist 
					`othvars_`demo''
					{ ;
			
					xi: reg `var' n_post `poly_`poly''
						`controls_`cov''
						if abs(dobdd_run)<=`bwidth' & keep_1997_2017_2c == 1 & demo_`demo' == 1
						, robust; 
					outreg2 using "`output'\\`exp'_bw`bwidth'_crime_poly`poly'`cov'cov_rd_`demo'.xls", e(N r2 rmse) append;
				} ; /* end var */		
				

		
		/* IV estimates */
				
				ivreg2 n_post (n_age18unfav=n_post) `poly_`poly''
					if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1
					, robust;
					outreg2 using  "`output'\\`exp'_bw`bwidth'_crime_poly`poly'`cov'cov_iv_`demo'.xls", e(N r2 rmse) replace;
				
				/* IV using unfavorable age 18 */
				foreach var of varlist 
					`vars_`demo''
					{;
					
					qui summ `var' 
						if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1 ;
					if r(sd) != 0 { ;
						xi: ivreg2 `var' (n_age18unfav=n_post) `poly_`poly''
							`controls_`cov''
							if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1
							, robust;
						outreg2 using  "`output'\\`exp'_bw`bwidth'_crime_poly`poly'`cov'cov_iv_`demo'.xls", e(N r2 rmse) append;
						} ;
					};
			
				if "`demo'" == "all" | "`demo'" == "male" | "`demo'" == "female" { ;
					foreach var of varlist 
						`hh_vars' 
						{ ;
					
						qui summ `var' 
							if abs(dobdd_run)<=`bwidth' & keep_1997_2017_2c == 1  & demo_`demo' == 1 ;
						if r(sd) != 0 { ;
							xi: ivreg2 `var' (n_age18unfav=n_post) `poly_`poly''
								`controls_`cov''
								if abs(dobdd_run)<=`bwidth' & keep_1997_2017_2c == 1  & demo_`demo' == 1
								, robust;
							outreg2 using  "`output'\\`exp'_bw`bwidth'_crime_poly`poly'`cov'cov_iv_`demo'.xls", e(N r2 rmse) append;
						} ;
					} ; /* end var */
				} ; /* end if */
				
			} ; /* end cov */
				
		} ; /* end poly */
	} ; /* end bwidth*/
} ;	/* end demo */



/* ******************************************** */
/* ***** YEAR-BY-YEAR ESTIMATES (STATSBY) ***** */
/* ******************************************** */
forval yr = 1997(1)2017 { ;
	rename yr_crime2genr_`yr' 	yr_genr_`yr' ;
	rename n_yrcrime2ahinc_`yr' yr_ahinc_`yr' ;
	rename n_yrcrime2shinc_`yr' yr_shinc_`yr' ;
	rename n_yrkidemp15_`yr' 	yr_kidemp15_`yr' ;
	rename n_yrkidearn`yr' 		yr_kidearn_`yr' ;
	rename n_onssi_`yr' 		yr_onssi_`yr' ;
} ;

local yr_list "genr ahinc shinc kidemp15 kidearn onssi" ;

foreach demo in all { ;
foreach bwidth in `bwidth_list' {;
forval poly = `poly_list' {;
foreach cov in all { ;

	foreach charge in `yr_list' { ;
		
		forval yr = 1997(1)2017 { ;

			xi: statsby _b _se, saving("`output'\statsby\\`exp'_yrbyyr_`charge'_`yr'.dta", replace):
				ivreg2 yr_`charge'_`yr' (n_age18unfav=n_post) `poly_`poly''
				`controls_`cov''
				if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1
				, robust;
						
			} ;
		
		} ;
} ;
} ;
} ;
} ;

preserve ;
foreach demo in all { ;
foreach bwidth in `bwidth_list' {;
forval poly = `poly_list' {;
foreach cov in all { ;
	
	foreach charge in `yr_list' { ;
		
		use "`output'\statsby\\`exp'_yrbyyr_`charge'_1997.dta", clear ;
		gen year = 1997 ;
		
		forval yr = 1997(1)2017 { ;
			append using "`output'\statsby\\`exp'_yrbyyr_`charge'_`yr'.dta" ;
			replace year = `yr' if year == . ;
		} ;
		drop if _n == 1 ;
		
		order year _b_n_age18unfav _se_n_age18unfav ;
		save 			"`output'\\`exp'_yrbyyr_`charge'_compiled.dta", replace ;
	
	} ;
				
} ;
} ;
} ;
} ;
restore ;


/* *********************************************** */
/* ***** CHARGE THRESHOLDS (WITH EMPLOYMENT) ***** */
/* *********************************************** */
local ct_list "" ;
foreach bwidth in `bwidth_list' {;				
foreach demo in all { ;
forval poly = `poly_list' {;
foreach cov in all { ;
		
	foreach var of varlist ct_* { ;

		qui summ `var'
			if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1 ;
		if r(sd) != 0 { ;
			local ct_list "`ct_list' `var'" ;
			xi: statsby _b _se, saving("`output'\statsby\\`exp'_`var'_`bwidth'_poly`poly'`cov'cov_iv_`demo'.dta", replace):
				 ivreg2 `var' (n_age18unfav=n_post) `poly_`poly''
				`controls_`cov''
				if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1
				, robust;
			} ;	
		
		} ;
} ;
} ;
} ;
} ;

preserve ;
foreach bwidth in `bwidth_list' {;				
foreach demo in all { ;
forval poly = `poly_list' {;
foreach cov in all { ;
				
			use 			"`output'\statsby\\`exp'_ct_crime2char1_`bwidth'_poly`poly'`cov'cov_iv_`demo'.dta", clear ;
			gen var = "ct_crime2char1" ;
		
			foreach var in `ct_list' { ;
				
				append using "`output'\statsby\\`exp'_`var'_`bwidth'_poly`poly'`cov'cov_iv_`demo'.dta" ;
				replace var = "`var'" if var == "" ;
					
				} ;
			} ;
			drop if _n == 1 ;
			
			save 			"`output'\\`exp'_ctcompiled_`bwidth'_poly`poly'`cov'cov_iv_`demo'.dta", replace ;
		
} ;
} ;
} ;
restore ;


/* ********************************************** */
/* ***** ROBUSTNESS TO DIFFERENT BANDWIDTHS ***** */
/* ********************************************** */
foreach var of varlist `robust_vars' {;
				
	foreach demo in all { ;

		forval poly = `poly_list' {;
			
			foreach cov in `cov_list' { ;
	
				forval bwidth = 25(25)400 { ;
				
					xi: statsby _b _se, saving("`output'\statsby\\`exp'_robust`bwidth'_`var'_poly`poly'`cov'cov_rd_`demo'.dta", replace):
						reg `var' n_post `poly_`poly''
							`controls_`cov''
							if abs(dobdd_run)<=`bwidth' & demo_`demo' == 1
							, robust; 
					} ;	
			
			} ;
		} ;
	} ;
} ;

preserve ;
foreach var of varlist `robust_vars' {;
				
	foreach demo in all { ;

		forval poly = `poly_list' {;
			
			foreach cov in `cov_list' { ;
			
				use 			"`output'\statsby\\`exp'_robust25_`var'_poly`poly'`cov'cov_rd_`demo'.dta", clear ;
				gen bwidth = 25 ;
				
				forval bwidth = 25(25)400 { ;
					append using "`output'\statsby\\`exp'_robust`bwidth'_`var'_poly`poly'`cov'cov_rd_`demo'.dta" ;
					replace bwidth = `bwidth' if bwidth == . ;
				} ;
				drop if _n == 1 ;
				
				order bwidth _b_n_post _se_n_post ;
				save 			"`output'\\`exp'_robust_`var'_poly`poly'`cov'cov_rd_`demo'.dta", replace ;
				
			} ;
		} ;
	} ;
} ;
restore ;

capture log close;
